利用{rayshader}生成全球任意位置的三维地形图. refer to https://github.com/camartinezbu/30DayChartChallenge2022/blob/f7f7cb1337f2604da2b78d7c920cb318606cf4c6/14-3-dimensional/scripts/1-plot.R
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ──────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6 ✔ purrr 0.3.4
✔ tibble 3.1.7 ✔ dplyr 1.0.9
✔ tidyr 1.2.0 ✔ stringr 1.4.0
✔ readr 2.1.2 ✔ forcats 0.5.1
── Conflicts ─────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
library(osmdata)
Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library(sf)
Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUE
library(rayshader)
library(raster)
载入需要的程辑包:sp
载入程辑包:‘raster’
The following object is masked from ‘package:dplyr’:
select
library(elevatr)
定义地理区域范围.
# map center longitude and latitude
center_lon_lat <- c(102.929, 30.150)
half_width <- 0.05
center_lon_lat <- c(116.321,39.947)
half_width <- 0.1
med_bbox <- st_bbox(c(xmin = center_lon_lat[1]-half_width, xmax = center_lon_lat[1]+half_width,
ymin = center_lon_lat[2]-half_width, ymax = center_lon_lat[2]+half_width),
crs = 4326)
med_bbox_df <- data.frame(x = c(center_lon_lat[1]-half_width, center_lon_lat[1]+half_width),
y = c(center_lon_lat[2]-half_width, center_lon_lat[2]+half_width))
extent_zoomed <- raster::extent(med_bbox)
利用{elevatr}下载地形数据.
prj_dd <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
elev_med <- get_elev_raster(med_bbox_df, prj=prj_dd, z=12, clip="bbox")
Accessing raster elevation [-------------------------] 0%
Accessing raster elevation [=>-----------------------] 8%
Accessing raster elevation [===>---------------------] 17%
Accessing raster elevation [=====>-------------------] 25%
Accessing raster elevation [=======>-----------------] 33%
Accessing raster elevation [=========>---------------] 42%
Accessing raster elevation [===========>-------------] 50%
Accessing raster elevation [==============>----------] 58%
Accessing raster elevation [================>--------] 67%
Accessing raster elevation [==================>------] 75%
Accessing raster elevation [====================>----] 83%
Accessing raster elevation [======================>--] 92%
Accessing raster elevation [=========================] 100%
Mosaicing & Projecting
Clipping DEM to bbox
Note: Elevation units are in meters.
elev_med_mat <- raster_to_matrix(elev_med)
[1] "Dimensions of matrix are: 1359x1358"
从OpenStreetMap获得地理信息数据
# Get highway data
# https://wiki.openstreetmap.org/wiki/Map_features#Highway
med_highway <- med_bbox %>% opq() %>%
add_osm_feature("highway",
c("motorway", "trunk", "primary", "secondary", "tertiary")) %>%
osmdata_sf()
med_highway_lines <- med_highway$osm_lines
# Get river data
# https://wiki.openstreetmap.org/wiki/Map_features#Waterway
med_rivers <- med_bbox %>% opq() %>%
add_osm_feature("waterway",
c("river", "stream")) %>%
osmdata_sf()
med_river_lines <- med_rivers$osm_lines
# Get buildings data
med_buildings <- med_bbox %>% opq() %>%
add_osm_feature("building") %>%
osmdata_sf()
med_building_polygons <- med_buildings$osm_polygons
# Get parkings data
med_parkdings <- med_bbox %>% opq() %>%
add_osm_feature("parking") %>%
osmdata_sf()
med_parking_polygons <- med_parkdings$osm_polygons
# Get tourism data
med_tourisms <- med_bbox %>% opq() %>%
add_osm_feature("tourism") %>%
osmdata_sf()
med_tourism_polygons <- med_tourisms$osm_polygons
# Get place data
med_places <- med_bbox %>% opq() %>%
add_osm_feature("place") %>%
osmdata_sf()
med_places_points <- med_places$osm_points %>%
filter(place %in% c("town","city"))
base_map <- elev_med_mat %>%
height_shade() %>%
add_overlay(sphere_shade(elev_med_mat, texture = "desert", colorintensity = 5), alphalayer=0.5) %>%
add_shadow(lamb_shade(elev_med_mat), 0) %>%
add_shadow(ambient_shade(elev_med_mat),0) %>%
add_shadow(texture_shade(elev_med_mat,detail=8/10,contrast=9,brightness = 11), 0.1)
plot_map(base_map)
# create highway layer
med_road_01 <- med_highway_lines %>%
filter(highway %in% c("motorway","trunk"))
med_road_02 <- med_highway_lines %>%
filter(highway %in% c("primary", "secondary"))
med_road_03 <- med_highway_lines %>%
filter(highway %in% c("tertiary"))
road_layer <- generate_line_overlay(med_road_03,extent = extent_zoomed,
linewidth = 2, color="white",
heightmap = elev_med_mat) %>%
add_overlay(generate_line_overlay(med_road_02,extent = extent_zoomed,
linewidth = 1, color="white", lty=3,
heightmap = elev_med_mat)) %>%
add_overlay(generate_line_overlay(med_road_01,extent = extent_zoomed,
linewidth = 3, color="white",
heightmap = elev_med_mat))
# create wateray layer
med_river_01 <- med_river_lines %>%
filter(waterway %in% c("river"))
med_river_02 <- med_river_lines %>%
filter(waterway %in% c("stream"))
water_layer <- generate_line_overlay(med_river_01,extent = extent_zoomed,
linewidth = 2, color="skyblue2",
heightmap = elev_med_mat) %>%
add_overlay(generate_line_overlay(med_river_02,extent = extent_zoomed,
linewidth = 1, color="skyblue2", lty=3,
heightmap = elev_med_mat))
# polygon_layer = generate_polygon_overlay(med_parking_polygons, extent = extent_zoomed,
# heightmap = elev_med_mat, palette="grey30") %>%
# add_overlay(generate_polygon_overlay(med_building_polygons, extent = extent_zoomed,
# heightmap = elev_med_mat, palette="darkred")) %>%
# add_overlay(generate_polygon_overlay(med_tourism_polygons, extent = extent_zoomed,
# heightmap = elev_med_mat, palette="darkgreen"), alphalayer = 0.6)
polygon_layer = generate_polygon_overlay(med_building_polygons, extent = extent_zoomed,
heightmap = elev_med_mat, palette="darkred")
Error in s2_geography_from_wkb(x, oriented = oriented, check = check) :
Evaluation error: Found 1 feature with invalid spherical geometry.
[28669] Loop 0 is not valid: Edge 7 crosses edge 10.
base_map %>%
# add_overlay(polygon_layer) %>%
add_overlay(water_layer, alphalayer = 0.8) %>%
add_overlay(road_layer) %>%
add_overlay(generate_label_overlay(med_places_points, extent = extent_zoomed,
text_size = 1, point_size = 1, color = "white",
halo_color = "black", halo_expand = 10,
halo_blur = 20, halo_alpha = 0.8,
seed=1, heightmap = elev_med_mat, data_label_column = "name.en")) %>%
plot_map()
base_map %>%
# add_overlay(polygon_layer) %>%
add_overlay(water_layer, alphalayer = 0.8) %>%
add_overlay(road_layer) %>%
add_overlay(generate_label_overlay(med_places_points, extent = extent_zoomed,
text_size = 1, point_size = 1, color = "white",
halo_color = "black", halo_expand = 10,
halo_blur = 15, halo_alpha = 0.8,
seed=1, heightmap = elev_med_mat, data_label_column = "name.en")) %>%
plot_3d(elev_med_mat, windowsize=c(1200,1100), zscale=20, theta=20, phi=30, fov=45, zoom=0.6)